﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Text.RegularExpressions;
using CalculateModelClasses;
using System.IO;
using TerFileProceed;


namespace RayTracingModel
{
   public class Intersection
    {
            //public 方法
            #region
                //求最近点时的比较函数

                public static int Compare(Pfd v1, Pfd v2)
                {
                    return v1.distance.CompareTo(v2.distance);
                }

                //[1]求射线方向上第一个与地形相交的点坐标
                //需要优化计算速度
                static public List<Pfd> IsIntersect(string TerInfo, Point one, SpectVector ot)
                {
                    List<Pfd> Crosspointandface = new List<Pfd>();
                    List<Pfd> crosspfd = new List<Pfd>();
                    List<KeyValuePair<string, List<double>>> material = TerMaterial.GetMaterial(TerInfo);
                    double a, b, c,conduct,permit;
                    string materialnum;
                    Point X, Y, Z;
                    try
                    {
                        
                        Regex r = new Regex("begin_<face>");
                        MatchCollection ma = r.Matches(TerInfo);
                        Regex p = new Regex("nVertices 3");
                        MatchCollection mc = p.Matches(TerInfo);
                        int i;
                        for (i = 0; i < mc.Count; i++)
                        {
                            int count = 0;
                            int m = mc[i].Index;
                            while (TerInfo[m + 13 + count] != ' ')
                            {
                                count++;
                            }
                            a = Convert.ToDouble(TerInfo.Substring(m + 13, count));
                            m = m + 13 + count + 1;
                            count = 0;
                            while (TerInfo[m + count] != ' ')
                            {
                                count++;
                            }
                            b = Convert.ToDouble(TerInfo.Substring(m, count));
                            m = m + count + 1;
                            count = 0;
                            while (TerInfo[m + count] != '\r')
                            {
                                count++;
                            }
                            c = Convert.ToDouble(TerInfo.Substring(m, count));
                            X = new Point(a, b, c);
                            m = m + count + 2;
                            count = 0;
                            while (TerInfo[m + count] != ' ')
                            {
                                count++;
                            }
                            a = Convert.ToDouble(TerInfo.Substring(m, count));
                            m = m + count + 1;
                            count = 0;
                            while (TerInfo[m + count] != ' ')
                            {
                                count++;
                            }
                            b = Convert.ToDouble(TerInfo.Substring(m, count));
                            m = m + count + 1;
                            count = 0;
                            while (TerInfo[m + count] != '\r')
                            {
                                count++;
                            }
                            c = Convert.ToDouble(TerInfo.Substring(m, count));
                            Y = new Point(a, b, c);
                            m = m + count + 2;
                            count = 0;
                            while (TerInfo[m + count] != ' ')
                            {
                                count++;
                            }
                            a = Convert.ToDouble(TerInfo.Substring(m, count));
                            m = m + count + 1;
                            count = 0;
                            while (TerInfo[m + count] != ' ')
                            {
                                count++;
                            }
                            b = Convert.ToDouble(TerInfo.Substring(m, count));
                            m = m + count + 1;
                            count = 0;
                            while (TerInfo[m + count] != '\r')
                            {
                                count++;
                            }
                            c = Convert.ToDouble(TerInfo.Substring(m, count));
                            Z = new Point(a, b, c);
                            int length = FileObject.Tool.readline(TerInfo, ma[i].Index);
                            int start = ma[i].Index + length + 2;
                            length = FileObject.Tool.readline(TerInfo, start);
                            materialnum = TerInfo.Substring(start, length);
                            int j;
                            for (j = 0; j < material.Count; j++)
                            {
                                if (material[j].Key == materialnum)
                                    break;
                            }
                            if (j >= material.Count) throw new Exception(" no respond material");
                            else
                            {
                                conduct = material[j].Value[0];
                                permit = material[j].Value[1];
                            }
                            Point CrossPoint = IsIntersectTri(one, ot, X, Y, Z);
                            Pfd temp;
                            Face correspond=new Face(X,Y,Z,conduct,permit);
                            if (CrossPoint != null)//如果返回point，则检测是否在三角面内
                            {
                                if (IsInthetri(X, Y, Z, CrossPoint))
                                {
                                    //if (Isthevertical(X, Y, Z, CrossPoint, e))
                                    //{
                                    //    Console.WriteLine("交点{0} {1} {2}在面{3}的顶点上", CrossPoint.x, CrossPoint.y, CrossPoint.z, i);
                                    //}
                                    //else if (OntheLine(X, Y, Z, CrossPoint, e, f))
                                    //{
                                    //    Console.WriteLine("交点{0} {1} {2}在面{3}的边上", CrossPoint.x, CrossPoint.y, CrossPoint.z, i);
                                    //}
                                    //else Console.WriteLine("交点{0} {1} {2}在面{3}的内部", CrossPoint.x, CrossPoint.y, CrossPoint.z, i);
                                    if ((CrossPoint.x - one.x) / ot.a > 0)
                                    {
                                        temp = new Pfd(CrossPoint,correspond, Distance(CrossPoint, one));
                                        //CrossPoints.Add(CrossPoint);
                                        Crosspointandface.Add(temp);
                                    }
                                    else continue;
                                }
                                else continue;
                            }
                            else continue;
                        }
                    }
                    catch (Exception e)
                    {
                        throw new Exception(e.Message + "地形出错");
                    }
                    if (Crosspointandface.Count != 0)
                    {
                        Crosspointandface.Sort(Compare);
                        if (Crosspointandface.Count >= 2)
                        {
                            if (Math.Abs(Crosspointandface[0].distance - Crosspointandface[1].distance )< 0.00000001)
                            {
                                crosspfd.Add(Crosspointandface[0]);
                                crosspfd.Add(Crosspointandface[1]);
                                return crosspfd;
                            }
                            else
                            {
                                crosspfd.Add(Crosspointandface[0]);
                                return crosspfd;
                            }
                        }
                        else
                        {
                            crosspfd.Add(Crosspointandface[0]);
                            return crosspfd;
                        }
                    }
                    else return crosspfd;
                }

                //[2]求两点的距离的函数

                public static double Distance(Point a, Point b)
                {
                    double distance;
                    distance = Math.Sqrt(Math.Pow((a.x - b.x), 2) + Math.Pow((a.y - b.y), 2) + Math.Pow((a.z - b.z), 2));
                    return distance;
                }

                //[3]求射线与面的交点的坐标

                public static Point IsIntersectTri(Point one, SpectVector ot, Point a, Point b, Point c)
                {
                    //先求三角面的法向量abXac
                    double TriD, temp1, temp2, t;
                    Point CrossPoint;
                    double[] ab = { b.x - a.x, b.y - a.y, b.z - a.z };
                    double[] ac = { c.x - a.x, c.y - a.y, c.z - a.z };
                    double[] abXac =   {// 
                        ab[1] * ac[2] - ab[2] * ac[1], // 
                        -ab[0] * ac[2] + ab[2] * ac[0],// 
                        ab[0] * ac[1] - ab[1] * ac[0],   };
                    TriD = -(abXac[0] * a.x + abXac[1] * a.y + abXac[2] * a.z);
                    temp1 = abXac[0] * one.x + abXac[1] * one.y + abXac[2] * one.z + TriD;
                    temp2 = abXac[0] * ot.a + abXac[1] * ot.b + abXac[2] * ot.c;
                    if (Math.Abs(temp2) < 0.00000001)
                    {
                        return null;
                    }
                    else
                    {
                        t = -(temp1 / temp2);
                        CrossPoint = new Point(ot.a * t + one.x, ot.b * t + one.y, ot.c * t + one.z);
                        return CrossPoint;
                    }
                }

                //[4]求射线与面的交点的坐标

                public static Point IsIntersectTri(Point one, SpectVector ot, Point a, SpectVector Vertical)
                {
                    double TriD, temp1, temp2, t;
                    Point CrossPoint;
                    //double[] ab = { b.x - a.x, b.y - a.y, b.z - a.z };
                    //double[] ac = { c.x - a.x, c.y - a.y, c.z - a.z };
                    //double[] abXac =   {// 
                    //        ab[1] * ac[2] - ab[2] * ac[1], // 
                    //        -ab[0] * ac[2] + ab[2] * ac[0],// 
                    //        ab[0] * ac[1] - ab[1] * ac[0],   };
                    TriD = -(Vertical.a * a.x + Vertical.b * a.y + Vertical.c * a.z);
                    temp1 = Vertical.a * one.x + Vertical.b * one.y + Vertical.c * one.z + TriD;
                    temp2 = Vertical.a * ot.a + Vertical.b * ot.b + Vertical.c * ot.c;
                    if (Math.Abs(temp2) < 0.00000001)
                    {
                        return null;
                    }
                    else
                    {
                        t = -(temp1 / temp2);
                        CrossPoint = new Point(ot.a * t + one.x, ot.b * t + one.y, ot.c * t + one.z);
                        return CrossPoint;
                    }
                }

                //[5]判断交点是否在三角面中函数

                public static bool IsInthetri(Point a, Point b, Point c, Point d)
                {
                    if (d == a || d == b || d == c)
                        return true;
                    else
                    {
                        //判断是否在三角面中
                        double[] ab = { b.x - a.x, b.y - a.y, b.z - a.z };
                        double[] ad = { d.x - a.x, d.y - a.y, d.z - a.z };
                        double[] bc = { c.x - b.x, c.y - b.y, c.z - b.z };
                        double[] bd = { d.x - b.x, d.y - b.y, d.z - b.z };
                        double[] cd = { d.x - c.x, d.y - c.y, d.z - c.z };
                        double[] ca = { a.x - c.x, a.y - c.y, a.z - c.z };
                        double[] adXab = { ad[1] * ab[2] - ad[2] * ab[1], -ad[0] * ab[2] + ad[2] * ab[0], ad[0] * ab[1] - ad[1] * ab[0], };
                        double[] bdXbc = { bd[1] * bc[2] - bd[2] * bc[1], -bd[0] * bc[2] + bd[2] * bc[0], bd[0] * bc[1] - bd[1] * bc[0], };
                        double[] cdXca = { cd[1] * ca[2] - cd[2] * ca[1], -cd[0] * ca[2] + cd[2] * ca[0], cd[0] * ca[1] - cd[1] * ca[0], };
                        double m = adXab[0] * bdXbc[0] + adXab[1] * bdXbc[1] + adXab[2] * bdXbc[2];
                        double n = bdXbc[0] * cdXca[0] + bdXbc[1] * cdXca[1] + bdXbc[2] * cdXca[2];
                        double l = cdXca[0] * adXab[0] + cdXca[1] * adXab[1] + cdXca[2] * adXab[2];
                        if (m >=0 && n >=0 && l >= 0)
                            return true;
                        else return false;
                    }
                }

                //[6]求三角面方程

                public static List<double> TriFoundation(Point a, Point b, Point c)
                {
                    List<double> para = new List<double>();
                    double D;
                    SpectVector paras = GetVerticalVector(a, b, c);
                    D = -(paras.a * a.x + paras.b * a.y + paras.c * a.z);
                    para.Add(paras.a); para.Add(paras.b); para.Add(paras.c);
                    para.Add(D);
                    return para;
                }

                //[7]求三角面方程

                public static List<double> TriFoundation(Face face)
                {
                    return TriFoundation(face.X, face.Y, face.Z);
                }

                //[8]判断点是否在面的外侧

                public static int IsOutFace(Face face, Point one)
                {
                    List<double> para = TriFoundation(face);
                    double temp = para[0] * one.x + para[1] * one.y + para[2] * one.z + para[3];
                    if (temp - 0 > 0.0001) return 2;
                    else if (temp - 0 < -0.0001) return 1;
                    else return 0;
                }

                //[9]判断点是否在面的外侧

                public static int IsOutFace(Point a, Point b, Point c, Point one)
                {
                    List<double> para = TriFoundation(a, b, c);
                    double temp = para[0] * one.x + para[1] * one.y + para[2] * one.z + para[3];
                    if (temp-0 > 0.0001) return 2;
                    else if (temp-0 < -0.0001) return 1;
                    else return 0;
                }

                //[10]求面的法向量

                public static SpectVector GetVerticalVector(Face face)
                {
                    return GetVerticalVector(face.X, face.Y, face.Z);
                }

                //[11]求面的法向量

                public static SpectVector GetVerticalVector(Point a, Point b, Point c)
                {
                    double[] ab = { b.x - a.x, b.y - a.y, b.z - a.z };
                    double[] ac = { c.x - a.x, c.y - a.y, c.z - a.z };
                    double[] abXac =   {// 
                        ab[1] * ac[2] - ab[2] * ac[1], // 
                        -ab[0] * ac[2] + ab[2] * ac[0],// 
                        ab[0] * ac[1] - ab[1] * ac[0],   };
                    SpectVector temp = new SpectVector(abXac[0], abXac[1], abXac[2]);
                    return temp;
                }

                //判断三个点是否在一条线段上的函数

                public static bool OnLine(Point Porta, Point Portb, Point c, double d)
                {
                    if (Math.Abs(Distance(Porta, c) + Distance(Portb, c) - Distance(Porta, Portb)) <= d)
                        return true;
                    else return false;
                }

                //判断该点是否在ter范围内

                public static bool Isinter(string TerInfo, Point one)
                {
                    int m = TerInfo.IndexOf("nVertices 3");
                    int count = 0;
                    while (TerInfo[m + 13 + count] != ' ')
                    {
                        count++;
                    }
                    double a = Convert.ToDouble(TerInfo.Substring(m + 13, count));
                    m = m + 13 + count + 1;
                    count = 0;
                    while (TerInfo[m + count] != ' ')
                    {
                        count++;
                    }
                    double b = Convert.ToDouble(TerInfo.Substring(m, count));
                    m = m + count + 1;
                    count = 0;
                    while (TerInfo[m + count] != '\r')
                    {
                        count++;
                    }
                    double c = Convert.ToDouble(TerInfo.Substring(m, count));
                    Point min = new Point(a, b, c);
                    m = TerInfo.LastIndexOf("nVertices 3") + 13;
                    count = 0;
                    while (TerInfo[m + count] != '\n')
                    {
                        count++;
                    }
                    m = m + count + 1;
                    count = 0;
                    while (TerInfo[m + count] != ' ')
                    {
                        count++;
                    }
                    a = Convert.ToDouble(TerInfo.Substring(m, count));
                    m = m + count + 1;
                    count = 0;
                    while (TerInfo[m + count] != ' ')
                    {
                        count++;
                    }
                    b = Convert.ToDouble(TerInfo.Substring(m, count));
                    m = m + count + 1;
                    count = 0;
                    while (TerInfo[m + count] != '\r')
                    {
                        count++;
                    }
                    c = Convert.ToDouble(TerInfo.Substring(m, count));
                    Point max = new Point(a, b, c);
                    if (one.x < min.x || one.x > max.x || one.y < min.y || one.y > max.y)
                        return false;
                    else return true;
                }

                //判断一个点是否在地形面上,返回所在面信息

                public static Face Isonter(string TerInfo, Point one)
                {
                    double a, b, c;
                    Point X, Y, Z;
                    Face temp = null;
                    List<Point> verticals = new List<Point>();
                    List<KeyValuePair<string, List<double>>> material = TerMaterial.GetMaterial(TerInfo);
                    string materialnum;
                    double conduct, permit;
                    Regex p = new Regex("begin_<face>");
                    MatchCollection ma = p.Matches(TerInfo);
                    Regex r = new Regex("nVertices 3");
                    MatchCollection mc = r.Matches(TerInfo);
                    int i;
                    for (i = 0; i < mc.Count; i++)
                    {
                        int count = 0;
                        int m = mc[i].Index;
                        while (TerInfo[m + 13 + count] != ' ')
                        {
                            count++;
                        }
                        a = Convert.ToDouble(TerInfo.Substring(m + 13, count));
                        m = m + 13 + count + 1;
                        count = 0;
                        while (TerInfo[m + count] != ' ')
                        {
                            count++;
                        }
                        b = Convert.ToDouble(TerInfo.Substring(m, count));
                        m = m + count + 1;
                        count = 0;
                        while (TerInfo[m + count] != '\r')
                        {
                            count++;
                        }
                        c = Convert.ToDouble(TerInfo.Substring(m, count));
                        X = new Point(a, b, c);
                        m = m + count + 2;
                        count = 0;
                        while (TerInfo[m + count] != ' ')
                        {
                            count++;
                        }
                        a = Convert.ToDouble(TerInfo.Substring(m, count));
                        m = m + count + 1;
                        count = 0;
                        while (TerInfo[m + count] != ' ')
                        {
                            count++;
                        }
                        b = Convert.ToDouble(TerInfo.Substring(m, count));
                        m = m + count + 1;
                        count = 0;
                        while (TerInfo[m + count] != '\r')
                        {
                            count++;
                        }
                        c = Convert.ToDouble(TerInfo.Substring(m, count));
                        Y = new Point(a, b, c);
                        m = m + count + 2;
                        count = 0;
                        while (TerInfo[m + count] != ' ')
                        {
                            count++;
                        }
                        a = Convert.ToDouble(TerInfo.Substring(m, count));
                        m = m + count + 1;
                        count = 0;
                        while (TerInfo[m + count] != ' ')
                        {
                            count++;
                        }
                        b = Convert.ToDouble(TerInfo.Substring(m, count));
                        m = m + count + 1;
                        count = 0;
                        while (TerInfo[m + count] != '\r')
                        {
                            count++;
                        }
                        c = Convert.ToDouble(TerInfo.Substring(m, count));
                        Z = new Point(a, b, c);
                        int length = FileObject.Tool.readline(TerInfo, ma[i].Index + 15);
                        materialnum = TerInfo.Substring(ma[i].Index + 15, length);
                        int j;
                        for (j = 0; j < material.Count; j++)
                        {
                            if (material[j].Key == materialnum)
                                break;
                        }
                        conduct = material[j].Value[0];
                        permit = material[j].Value[1];
                        temp = new Face(X, Y, Z, conduct, permit);
                        if (one.x < X.x || one.x > Y.x || one.y < X.y || one.y > Z.y)//判断是否在该三角面范围中
                        {
                            temp = null;
                            continue;
                        }
                        else
                        {
                            //if (isInOnePlane(X, Y, Z, one)) //如果四个点在一个平面上，判断该点是否在该三角面中
                            //{
                            //    if (isInthetri(X, Y, Z, one))
                            //    {
                            //        Console.WriteLine("在地形面上");
                            //        break;
                            //    }
                            //    else continue;
                            //}
                            //else continue;
                            //Zvalue(X, Y, Z, one);
                            if (!IsInthetri(X, Y, Z, one))
                            {
                                temp = null;
                                continue;
                            }
                        }

                    }
                    return temp;

                }

                //找到三角面的非棱点

                public static Point Find(Point Porta, Point Portb, Face face)
                {
                    if (face.X.unequal(Porta) && face.X.unequal(Portb)) return face.X;
                    else if (face.Y.unequal(Porta) && face.Y.unequal(Portb)) return face.Y;
                    else return face.Z;
                }

                //判断射线始点在两个相交面的位置，以及判断两个面的形状

                public static bool Iscontex(Point Porta, Point Portb, Point one, List<Face> faces)
                {
                    Point p0 = Find(Porta, Portb, faces[0]); Point p1 = Find(Porta, Portb, faces[1]);
                    int p0top1 = Intersection.IsOutFace(faces[1], p0); int p1top0 = Intersection.IsOutFace(faces[0], p1);
                    int txtoface0 = Intersection.IsOutFace(faces[0], one); int txtoface1 = Intersection.IsOutFace(faces[1], one);
                    if (txtoface0 == p1top0 && txtoface1 == p0top1) return false;//说明两个面相对于点是凹型的
                    else return true;

                }

                public static void GetZValue(string TerInfo, Point one)
                {
                    double a, b, c;
                    Point X, Y, Z;
                    Regex p = new Regex("nVertices 3");
                    MatchCollection mc = p.Matches(TerInfo);
                    int i;
                    for (i = 0; i < mc.Count; i++)
                    {
                        int count = 0;
                        int m = mc[i].Index;
                        while (TerInfo[m + 13 + count] != ' ')
                        {
                            count++;
                        }
                        a = Convert.ToDouble(TerInfo.Substring(m + 13, count));
                        m = m + 13 + count + 1;
                        count = 0;
                        while (TerInfo[m + count] != ' ')
                        {
                            count++;
                        }
                        b = Convert.ToDouble(TerInfo.Substring(m, count));
                        m = m + count + 1;
                        count = 0;
                        while (TerInfo[m + count] != '\r')
                        {
                            count++;
                        }
                        c = Convert.ToDouble(TerInfo.Substring(m, count));
                        X = new Point(a, b, c);
                        m = m + count + 2;
                        count = 0;
                        while (TerInfo[m + count] != ' ')
                        {
                            count++;
                        }
                        a = Convert.ToDouble(TerInfo.Substring(m, count));
                        m = m + count + 1;
                        count = 0;
                        while (TerInfo[m + count] != ' ')
                        {
                            count++;
                        }
                        b = Convert.ToDouble(TerInfo.Substring(m, count));
                        m = m + count + 1;
                        count = 0;
                        while (TerInfo[m + count] != '\r')
                        {
                            count++;
                        }
                        c = Convert.ToDouble(TerInfo.Substring(m, count));
                        Y = new Point(a, b, c);
                        m = m + count + 2;
                        count = 0;
                        while (TerInfo[m + count] != ' ')
                        {
                            count++;
                        }
                        a = Convert.ToDouble(TerInfo.Substring(m, count));
                        m = m + count + 1;
                        count = 0;
                        while (TerInfo[m + count] != ' ')
                        {
                            count++;
                        }
                        b = Convert.ToDouble(TerInfo.Substring(m, count));
                        m = m + count + 1;
                        count = 0;
                        while (TerInfo[m + count] != '\r')
                        {
                            count++;
                        }
                        c = Convert.ToDouble(TerInfo.Substring(m, count));
                        Z = new Point(a, b, c);
                        if (Min(X.x, Y.x, Z.x) < one.x && one.x < Max(X.x, Y.x, Z.x) && Min(X.y, Y.y, Z.y) < one.y && one.y < Max(X.y, Y.y, Z.y))
                        {
                            Zvalue(X, Y, Z, one);
                            if (IsInthetri(X, Y, Z, one))
                                break;
                            else continue;
                        }
                        else continue;
                    }
                }
                #endregion

            //私有方法
            #region

            
            //判断交点是否在三角形边上的函数
            static bool OnTheLine(Point a, Point b, Point c, Point d, double e, double f)
            {
                if (OnLine(a, b, d, f) || OnLine(a, c, d, f) || OnLine(b, c, d, f))
                    return true;
                else return false;

            }
            //已知点的XY值和所在面求Z值的函数
            static void Zvalue(Point a, Point b, Point c, Point one)
            {
                double[] ab = { b.x - a.x, b.y - a.y, b.z - a.z };
                double[] ac = { c.x - a.x, c.y - a.y, c.z - a.z };
                double[] abXac =   {// 
                    ab[1] * ac[2] - ab[2] * ac[1], // 
                    -ab[0] * ac[2] + ab[2] * ac[0],// 
                    ab[0] * ac[1] - ab[1] * ac[0],   };
                one.z = (-(abXac[0] * (one.x - a.x) + abXac[1] * (one.y - a.y)) / abXac[2]) + a.z;

            }
            //求Z最大值
            static double Max(double a, double b, double c)
            {
                return Math.Max(Math.Max(a, b), Math.Max(b, c));
            }
            //求Z最小值   
            static double Min(double a, double b, double c)
            {
                return Math.Min(Math.Min(a, b), Math.Min(b, c));
            }
            //判断四个点是否在一个面上(abxac)*ad=0
            static bool IsInOnePlane(Point a, Point b, Point c, Point d)  
            {
                double[] ab = { b.x - a.x, b.y - a.y, b.z - a.z };
                double[] ac = { c.x - a.x, c.y - a.y, c.z - a.z };
                double[] ad = { d.x - a.x, d.y - a.y, d.z - a.z };
                double[] abXac =   {// 
                    ab[1] * ac[2] - ab[2] * ac[1], // 
                    -ab[0] * ac[2] + ab[2] * ac[0],// 
                    ab[0] * ac[1] - ab[1] * ac[0],   };
                return Math.Abs(abXac[0] * ad[0] + abXac[1] * ad[1] + abXac[2] * ad[2]) < 0.000001;
            }
            //判断交点是否为顶点的函数 
            static bool IsTheVertical(Point a, Point b, Point c, Point d, double e)
            {
                if (Distance(a, d) <= e || Distance(b, d) <= e || Distance(c, d) <= e)
                    return true;
                else return false;
            }
            #endregion

    }
}
